Mapping data in R

Setup

Welcome to another tutorial for this class, COMP/STAT 112: Introduction to Data Science! It will be similar to the others, including demo videos and files embedded in this document and practice problems with hints or solutions at the end. There are some new libraries, so be sure to install those first.

As most of our files do, we start this one with three R code chunks: 1. options, 2. libraries and settings, 3. data.

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

library(tidyverse)     # for data cleaning and plotting
library(googlesheets4) # for reading googlesheet data
library(lubridate)     # for date manipulation
library(openintro)     # for the abbr2state() function
library(palmerpenguins)# for Palmer penguin data
library(maps)          # for map data
library(ggmap)         # for mapping points on maps
library(gplots)        # for col2hex() function
library(RColorBrewer)  # for color palettes
library(sf)            # for working with spatial data
library(leaflet)       # for highly customizable mapping
library(ggthemes)      # for more themes (including theme_map())
gs4_deauth()           # To not have to authorize each time you knit.
theme_set(theme_minimal())

# Starbucks locations
Starbucks <- read_csv("https://www.macalester.edu/~ajohns24/Data/Starbucks.csv")

# Lisa's favorite St. Paul places - used in leaflet example
favorite_stp_by_lisa <- tibble(
  place = c("Home", "Macalester College", "Adams Spanish Immersion", 
            "Spirit Gymnastics", "Bama & Bapa", "Now Bikes",
            "Dance Spectrum", "Pizza Luce", "Brunson's"),
  long = c(-93.1405743, -93.1712321, -93.1451796, 
           -93.1650563, -93.1542883, -93.1696608, 
           -93.1393172, -93.1524256, -93.0753863),
  lat = c(44.950576, 44.9378965, 44.9237914,
          44.9654609, 44.9295072, 44.9436813, 
          44.9399922, 44.9468848, 44.9700727)
  )

Learning Goals

After this tutorial, you should be able to do the following:

• Plot data points on top of a map using the ggmap() function along with ggplot2 functions.

• Create choropleth maps using geom_map().

• Add points and other ggplot2 features to a map created from geom_map().

• Understand the basics of creating a map using leaflet, including adding points and choropleths to a base map.

Plotting points on a map

The Starbucks data, compiled by Danny Kaplan and provided by Alicia Johnson, contains information about every Starbucks in the world at the time the data were collected. It includes the Latitude and Longitude of each location. Let’s start by using familiar plotting tools

ggplot(data=Starbucks) +
  geom_point(aes(x = Longitude, y = Latitude), 
             alpha = 0.2, 
             size = .1)

The point pattern probably looks familiar. To highlight the geographical nature of this scatterplot, we can superimpose the points on top of a map, using the ggmap() function from the ggmap library.

NOTE: we used to be able to easily bring in Google maps. As of mid-2018, in order to bring those in, you need to have a registered API key. If you want to do that, see google_key in the help. Then, see the documentation for get_map(). We will bring in other types of maps since Google maps are harder to do now and require you to submit credit card information.

Instead, we bring in a stamen map (there are others you could try, but we’ll stick with this). You can also take a look at stamen maps on their website. First, let’s look at an example.

# Get the map information
world <- get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

# Plot the points on the map
ggmap(world) + # creates the map "background"
  geom_point(data = Starbucks, 
             aes(x = Longitude, y = Latitude), 
             alpha = .3, 
             size = .1) +
  theme_map()

Next, we will walk through the get_stamenmap() function arguments. The code below is what was used to get the world map information.

get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

bbox

get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

The bbox argument tells it the minimum and maximum latitude and longitude points. So, left is the minimum longitude, right is the maximum longitude, bottom is the minimum latitude, and top is the maximum latitude. I found it helpful to go to openstreetmap: zoom in on the area of interest, click export, and you will see all the values you need. I had to modify them slightly, which you can do after your initial plot.

maptype

get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

The maptype tells it the style of the map. Check out the different options by looking in the get_stamenmap help (type ?get_stamenmap in the console).

zoom

get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

When you make a large area, you need to decrease the zoom, otherwise it will take too long to load. So, it’s a good idea to start with a small zoom and you can always make it bigger if you want. This might seem counter-intuitive at first. I think of the zoom level as the level of detail. So, smaller numbers show less detail and larger numbers more detail. I often go to the stamanmaps webpage and search for the location I’m mapping. Then, in the URL, you can see the zoom number. For example, this link is a map of St. Paul: http://maps.stamen.com/#terrain/12/44.9531/-93.0904. Notice the number 12 next to /#terrain/. That means it is zoomed in at 12.

ggmap()

We save the the map information from get_stamenmap() to a named value and then use it in ggmap():

# Get the map information
world <- get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

# Plot the points on the map
ggmap(world) + # creates the map "background"
  geom_point(data = Starbucks, 
             aes(x = Longitude, y = Latitude), 
             alpha = .3, 
             size = .1) +
  theme_map()

The ggmap() function will print the “background” map. Think of it as the providing the canvas on which we will plot. This takes the place of our usual ggplot().

ggmap(world)

After that, we can use the geom_XXX() functions from ggplot2 that we are used to in order to put points, lines, etc. on top of the map. But, we need to remember to also provide the data we are using in the geom_XXX() function(s) we use since we do not have the ggplot() function in which to provide it.

# Get the map information
world <- get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

# Plot the points on the map
ggmap(world) + # creates the map "background"
  geom_point(data = Starbucks, 
             aes(x = Longitude, y = Latitude), 
             alpha = .3, 
             size = .1) +
  theme_map()

theme_map()

The last thing I did in the code was to add theme_map(). This is optional, but I often find it makes it look nice.

# Get the map information
world <- get_stamenmap(
    bbox = c(left = -180, bottom = -57, right = 179, top = 82.1), 
    maptype = "terrain",
    zoom = 2)

# Plot the points on the map
ggmap(world) + # creates the map "background"
  geom_point(data = Starbucks, 
             aes(x = Longitude, y = Latitude), 
             alpha = .3, 
             size = .1) +
  theme_map()

So, the final map as a world map as the background with points plotted on top that show the Starbucks locations. The points are .1 of their usual size and have a transparency level of .3.

Demo video

With that introduction, you are ready to watch the demo video!

Voicethread: ggmap demo

Download ggmap demo file (without code) Download ggmap demo file (with code)

Resources

• Examples from ggmap maintainer David Kahle
  
• ggmap cheatsheet

Your turn!

Exercise: More with Starbucks

1. Add an aesthetic to the world map that sets the color of the points according to the ownership type. What, if anything, can you deduce from this visualization?

2. Construct a new map of Starbucks locations in the Twin Cities metro area (approximately the 5 county metro area).

3. In the Twin Cities plot, play with the zoom number. What does it do? (just describe what it does - don’t actually include more than one map).

4. Try a couple different map types (see get_stamenmap() in help and look at maptype). Include a map with one of the other map types.

5. Add a point to the map that indicates Macalester College and label it appropriately. There are many ways you can do think, but I think it’s easiest with the annotate() function (see ggplot2 cheatsheet).

Choropleths

Geographical data needn’t be expressed by latitude and longitude. For choropleth maps, instead of visualizing our data as points with different aesthetics (size, color, transparency, etc.), we color different regions of the maps based on data values. To do this we need to specify both the geometric regions on which the data resides (counties, states, zip codes, etc.), and then wrangle the data so that there is one value per region.

Let’s return to the Starbucks data. First, we will create a new dataset, starbucks_us_by_state that limits the data to the US, finds the number of Starbucks in each state, and creates a state name that is in all lowercase letters that matches the state name in the region variable of the states_map dataset.

The states_map dataset gives information about creating the borders of the US states. The data is retrieved using the map_data() function. Run ?map_data in the console to see more information about what other maps are available. There are also other packages that provide different types of maps.

Then, we can use geom_map() to create a choropleth map. Let’s take a look at the map and we’ll go through the details after.

#Create a new Starbucks dataset that 
# - filters to the US
# - summarizes the number of Starbucks in each state
# - has full names of states in lowercase letters (to match to states_map data created next)

starbucks_us_by_state <- Starbucks %>% 
  filter(Country == "US") %>% 
  count(`State/Province`) %>% 
  mutate(state_name = str_to_lower(abbr2state(`State/Province`))) 

#US states map information - coordinates used to draw borders
states_map <- map_data("state")

# map that colors state by number of Starbucks
starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  #This assures the map looks decently nice:
  expand_limits(x = states_map$long, y = states_map$lat) + 
  theme_map()

Now, let’s look more closely at what each piece of the code below is doing.

starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  expand_limits(x = states_map$long, y = states_map$lat) + 
  theme_map()

Choose a map

The map argument tells R at which level to create the map. Really, it tells it how to draw all the borders This is a very special data set. According to the geom_map() documentation, it is a “data frame that contains the map coordinates … It must contain columns x or long, y or lat, and region or id.” We are using the map_data() function to create the map file (see above for more detail0. You can open the map data, states_map, and see that it adheres to the rules.

starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  expand_limits(x = states_map$long, y = states_map$lat) + 
  theme_map()

Connect map id/region variable to data being plotted

The map_id inside of aes() is a required aesthetic for the geom_map() geom. It tells R which variable is the region/id variable, in this case the state. It connects the region or id from the map (region variable in states_map dataset, in this example) to the dataset being plotted (state_name in starbucks_us_by_state, in this example). So state_name needs to have the same form as region, which is why we modified the state names in starbucks_us_by_state.

starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  expand_limits(x = states_map$long, y = states_map$lat) + 
  theme_map()

Use ggplot2 features

We tell it to fill in the states by the variable n, the number of Starbucks in each state. With the geom_map() geom, it will fill in the borders of the regions we defined in the map argument.

starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  expand_limits(x = states_map$long, y = states_map$lat) + 
  theme_map()

expand_limits()

Use expand_limits() to assure that the map covers the entire area it’s supposed to. We put the longitude variable from states_map for the x argument and the latitude variable from states_map for the y argument to assure the map stretches across the entire range of longitudes and latitudes in the map. There may be a better way to do this but I have yet to find it, and when I leave it out, I don’t even see the map.

starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  expand_limits(x = states_map$long, y = states_map$lat) + 
  theme_map()

theme_map()

This is a personal preference. I like the way theme_map() makes the map look.

starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  expand_limits(x = states_map$long, y = states_map$lat) + 
  theme_map()

Add ggplot2 layers

You can add any of the ggplot2 layers on top of this map. In this example, I’ve added the MN Starbucks as points, a title, and changed the legend background (so it doesn’t have one and overlap California).

starbucks_us_by_state %>% 
  ggplot() +
  geom_map(map = states_map,
           aes(map_id = state_name,
               fill = n)) +
  geom_point(data = Starbucks %>% filter(`State/Province` == "MN"),
             aes(x = Longitude, y = Latitude),
             size = .05,
             alpha = .2, 
             color = "goldenrod") +
  expand_limits(x = states_map$long, y = states_map$lat) + 
  labs(title = "Starbucks in MN") +
  theme_map() +
  theme(legend.background = element_blank())

Demo video

And you are ready to watch another demo video!

Voicethread: geom_map() demo

Download geom_map() demo file (without code) Download geom_map() demo file (with code)

Resources

• ggplot2 documentation
  
• Example by Arie Voorman (some things could be out of date since it’s from 2015)

Your turn!

Exercise: Even more with Starbucks

The example I showed did not account for population of each state in the map. In the code below, a new variable is created, starbucks_per_10000, that gives the number of Starbucks per 10,000 people. It is in the starbucks_with_2018_pop_est dataset.

census_pop_est_2018 <- read_csv("https://www.dropbox.com/s/6txwv3b4ng7pepe/us_census_2018_state_pop_est.csv?dl=1") %>% 
  separate(state, into = c("dot","state"), extra = "merge") %>% 
  select(-dot) %>% 
  mutate(state = str_to_lower(state))

starbucks_with_2018_pop_est <-
  starbucks_us_by_state %>% 
  left_join(census_pop_est_2018,
            by = c("state_name" = "state")) %>% 
  mutate(starbucks_per_10000 = (n/est_pop_2018)*10000)

1. dplyr review: Look through the code above and describe what each line of code does.

2. Create a choropleth map that shows the number of Starbucks per 10,000 people on a map of the US. Use a new fill color, add points for all Starbucks in the US (except Hawaii and Alaska), add an informative title for the plot, and include a caption that says who created the plot (you!). Make a conclusion about what you observe.

Using leaflet to create maps

Concept Map

Leaflet is an open-source JavaScript library for creating maps. It can be used outside of R, but we will only discuss using the leaflet library in R.

This library uses a different plotting framework from ggplot2 although it still has a tidyverse feel due to its use of the pipe, %>% and the way it adds layers to the plot, just like in ggplot2.

Introductory video

Watch the video that introduces leaflet() functions.

Voicethread: Mapping in R with leaflet

Download slides

Steps to create a map

1. Create a map widget by calling leaflet() and telling it the data to use.
   
2. Add a base map using addTiles() (the default) or addProviderTiles().
3. Add layers to the map by using layer functions (e.g. , addMarkers(), addPolygons()) to modify the map widget.
   
4. Repeat step 3 as desired.
   
5. Print the map widget to display it.

Creating a map

Now, I create a basic map and add my points (the points are a layer on the map). The data are in favorite_stp_by_lisa.

The function we will use to create the maps will look for certain variable names for latitude (lat, latitude) and longitude (lng, long, or longitude). If you do not name them one of those things or if the data you are using doesn’t name them that, you need to call out the name explicitly. You can use a “two-finger scroll” to zoom in and out.

leaflet(data = favorite_stp_by_lisa) %>% #base plot
  addTiles() %>% #base map - default is openstreet map but that can be changed
  addMarkers() #Adds markers - knows lat and long from names in data

Same as above but explicitly told it latitude and longitude, which you would need to do if those variables had a name not recognized by the function, and added labels. WARNING: DO NOT FORGET THE ~ BEFORE THE VARIABLE NAMES!!!

leaflet(data = favorite_stp_by_lisa) %>% 
  addTiles() %>% 
  addMarkers(lng = ~long, lat = ~lat, label = ~place) 

We can change just about everything about our map. This is the same plot as above with some aesthetic changes. Some tips:

• To see all available provider base maps, type providers in the console.

• To access those maps, use providers$PROVIDERNAME inside the addProviderTiles() function, where PROVIDERNAME is one of those listed providers. When you type provider$ a list should show up that you can click on.

• Colors need to be in “hex” form. I used the col2hex() function from the gplot library to do that since I don’t have any hex colors memorized. You can also google hex colors and find websites to help you out.

• Search addControl in the Help or type ?addControl into the console to see what all the arguments mean and how you can change them.

leaflet(data = favorite_stp_by_lisa) %>% 
  addProviderTiles(providers$Stamen.Watercolor) %>% 
  addCircles(lng = ~long, 
             lat = ~lat, 
             label = ~place, 
             weight = 10, 
             opacity = 1, color = col2hex("darkblue")) 

The map below is also the “same” as the ones I have already created with a new base map and a line to trace my route, which was created with the addPolylines() layer. It traces it in the order they are entered in the dataset.

leaflet(data = favorite_stp_by_lisa) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addCircles(lng = ~long, 
             lat = ~lat, 
             label = ~place, 
             weight = 10, 
             opacity = 1, 
             color = col2hex("darkred")) %>% 
  addPolylines(lng = ~long, 
               lat = ~lat, 
               color = col2hex("darkred"))

Choropleth layers with addPolygons()

Resources

• Detailed documenation (with examples)

• Cheatsheet

• Provider map previews